# Equivalence test
# according to page 3 in https://cran.r-project.org/web/packages/equivalenceTest/equivalenceTest.pdf

my1<-shared_hits$beta_elisa
my2<-shared_hits$beta_olink_corrected
se1<-shared_hits$StdErr_elisa
se2<-shared_hits$StdErr_olink
#different Deltas can be tested 
delta<-0.333
delta<- 0.1


shared_hits$tau1<-(my2-my1+delta)/sqrt(se1^2+se2^2)
shared_hits$tau2<-(my2-my1-delta)/sqrt(se1^2+se2^2)
shared_hits$tau1bigger<-shared_hits$tau1>qnorm(0.95)
shared_hits$tau2smaller<-shared_hits$tau2<(-qnorm(0.95))

table(shared_hits$tau1bigger) 
table(shared_hits$tau2smaller) 


shared_hits[shared_hits$tau1bigger==F,]
shared_hits[shared_hits$tau2smaller==F,]



# adjusting for multiple testing 
shared_hits$tau1bigger<-shared_hits$tau1>qnorm(1-0.05/2110)
shared_hits$tau2smaller<-shared_hits$tau2<(-qnorm(1-0.05/2110))

table(shared_hits$tau1bigger) 
table(shared_hits$tau2smaller) 


shared_hits[shared_hits$tau1bigger==F,]  
shared_hits[shared_hits$tau2smaller==F,] 

shared_hits$lowerCI.elisa<-shared_hits$beta_elisa-qnorm(0.975)*shared_hits$StdErr_elisa
shared_hits$upperCI.elisa<-shared_hits$beta_elisa+qnorm(0.975)*shared_hits$StdErr_elisa

shared_hits$lowerCI.olink<-shared_hits$beta_olink_corrected-qnorm(0.975)*shared_hits$StdErr_olink
shared_hits$upperCI.olink<-shared_hits$beta_olink_corrected+qnorm(0.975)*shared_hits$StdErr_olink



#########################
cor.test(shared_hits$beta_olink_corrected,shared_hits$beta_elisa)

shared_hits$same_direction <- shared_hits$beta_elisa * shared_hits$beta_olink_corrected > 0

# Calculate SNPs with the same effect direction
shared_hits <- shared_hits %>%
  mutate(same_direction = beta_elisa * beta_olink_corrected > 0)

# Count total SNPs
total_snps <- nrow(shared_hits)

# Count SNPs where ELISA and Olink have the same effect direction
same_direction_snps <- sum(shared_hits$same_direction, na.rm = TRUE)

# Calculate percentage
same_direction_percentage <- (same_direction_snps / total_snps) * 100

# Print results
cat("Total Shared SNPs:", total_snps, "\n") 
cat("SNPs with Same Effect Direction:", same_direction_snps, "\n") 
cat("Percentage of SNPs with Same Effect Direction:", round(same_direction_percentage, 2), "%\n") 


################
#final plot

plot(shared_hits$beta_elisa,
     shared_hits$beta_olink_corrected,
     pch=19,
     col="darkgrey",
     xlim=c(-0.4,0.4),
     ylim=c(-0.4,0.4),
     xlab="",
     ylab="",
     bty="n")

# Add bold axis labels
title(xlab = "beta estimate +/- 95% CI from ELISA", font.lab = 2, cex.lab = 1.2)
title(ylab = "beta estimate +/- 95% CI from UKB", font.lab = 2, cex.lab = 1.2)

#labels for specific variants can be added manually to avoid overlap
# example Variant 1 – label above
points(shared_hits$beta_elisa[shared_hits$MarkerName == "7:81970590"],
       shared_hits$beta_olink_corrected[shared_hits$MarkerName == "7:81970590"],
       pch = 19, col = "darkblue")
text(shared_hits$beta_elisa[shared_hits$MarkerName == "7:81970590"]+ 0.125,
     shared_hits$beta_olink_corrected[shared_hits$MarkerName == "7:81970590"] - 0.025,
     labels = "rs79741398 (7:81970590) ", col = "darkblue", cex = 0.8)




abline(0,1,lty=2)
abline(lm(shared_hits$beta_elisa~shared_hits$beta_olink_corrected),lty=3)
abline(h=0,lty=1)
abline(v=0,lty=1)

segments(shared_hits$lowerCI.elisa,shared_hits$beta_olink_corrected,shared_hits$upperCI.elisa,shared_hits$beta_olink_corrected,col="darkgrey")
segments(shared_hits$beta_elisa,shared_hits$lowerCI.olink,shared_hits$beta_elisa,shared_hits$upperCI.olink,col="darkgrey")



highlighted_idx <- shared_hits$MarkerName %in% c("7:81971189", "13:55097589", "11:117111368", "11:116692293", 
                                                                 "11:117042408", "7:81958748", "7:81970590")

# Horizontal bars (x = CI of ELISA, y = Effect_olink)
segments(shared_hits$lowerCI.elisa[highlighted_idx],
         shared_hits$beta_olink_corrected[highlighted_idx],
         shared_hits$upperCI.elisa[highlighted_idx],
         shared_hits$beta_olink_corrected[highlighted_idx],
         col = "darkblue")

# Vertical bars (x = Effect_elisa, y = CI of Olink)
segments(shared_hits$beta_elisa[highlighted_idx],
         shared_hits$lowerCI.olink[highlighted_idx],
         shared_hits$beta_elisa[highlighted_idx],
         shared_hits$upperCI.olink[highlighted_idx],
         col = "darkblue")


shared_hits$highlight_blue <-
  !(shared_hits$tau1bigger & shared_hits$tau2smaller)


#exclude those that are already highlighted in darkblue (when margin 0.33 is used)
shared_hits$highlight_blue[highlighted_idx] <- FALSE
summary(shared_hits$highlight_blue)


with(subset(shared_hits, highlight_blue),
     points(beta_elisa, beta_olink_corrected, pch = 19, col = "gray50"))

